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We perform highly accurate density matrix renormalization group (DMRG) simulations to investigate the 
ground state properties of the spin-j antiferromagnetic square lattice Heisenberg J1-J2 model. Based on stud- 
ies of numerous long cylinders with circumferences of up to 14 lattice spacings, we obtain strong evidence for 
a topological quantum spin liquid state in the region 0.41 < J2/J1 < 0.62, separating conventional Neel and 
striped antiferromagnetic states for smaller and larger J%j J\, respectively. The quantum spin liquid is char- 
acterized numerically by the absence of magnetic or valence bond solid order, and non-zero singlet and triplet 
energy gaps. Furthermore, we positively identify its topological nature by measuring a non-zero topological 
entanglement entropy 7 = 0.70 ± 0.02, extremely close to 7 = ln(2) « 0.69 (expected for a Zi quantum 
spin liquid) and a non-trivial finite size dimerization effect depending upon the parity of the circumference of 
the cylinder. We also point out that a valence bond solid, and indeed any discrete symmetry breaking state, 
would be expected to show a constant correction to the entanglement entropy of opposite sign to the topological 
entanglement entropy. 

PACS numbers: 75.10.Jm, 75.50.Ee, 75.40.Mg 



I. INTRODUCTION 

Quantum spin liquids (QSLs) are elusive magnets without 
magnetism, resisting symmetry breaking even at zero tem- 
perature due to strong quantum fluctuations and geometric 
frustration 1 . The simplest QSLs known theoretically are char- 
acterized by topological order 2 ^, and support fractionalized 
excitations including spinons, which carry the spin (1/2) but 
not the charge of the electron. Since the QSL state was 
suggested by Anderson 5 , it has been sought, mostly unsuc- 
cessfully, in models and materials. However, exciting indi- 
cations of QSL ground states were recently reported in nu- 
merical studies of models on the honeycomb 6 and kagome 7 
lattices. Here we report strong evidence for a QSL state in 
the square lattice J\-Ji antiferromagnetic (AFM) Heisenberg 
model, with the Hamiltonian 
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where S, is the spin- 1/2 operator on site i and (ij) (((ij))) 
denotes nearest neighbors (next nearest neighbors). In the fol- 
lowing we set Ji — 1 as the unit of energy, and consider only 
the frustrated case J 2 > 0. 

Eq. ([TJ is of fundamental interest for its simplicity, and for 
its relevance to cuprates, Fe-based superconductor*^ 2 -! an( j 
other materials^. Accordingly, it is among the most studied 
models in frustrated quantum magnetism^ 2 !! These previ- 
ous studies have established the existence of a non-magnetic 
ground state between the Neel and striped AFM states which 
occur for small and large J 2, respectively. 



To characterize the non-magnetic phase, we can ask two 
main types of questions. First, we may ask about its sym- 
metries. Being non-magnetic, the ground state retains the in- 
ternal SU(2) spin-rotation invariance, but it may break spatial 
ones. If SU(2) is preserved but spatial symmetries are bro- 
ken in such a way that the unit cell is enlarged, the system 
is said to have valence bond solid (VBS) order. Second, we 
may ask about the range of entanglement of the wavefunction. 
The simplest representative wavefunctions for VBS states are 
continuously deformable by local unitary transformations into 
product states. Such is true for typical ground state wavefunc- 
tions for systems with broken discrete symmetries (the space 
group of a lattice is discrete). As such, these wavefunctions 
have only short-range entanglement (Schrodinger cat states 



are possible in finite systems and will be discussed in Sec. IV 1 



Wavefunctions which cannot be continuously transformed in 
this way into product states may be said to exhibit long-range 
entanglement. This is true for all gapless critical phases, as 
well as for some gapped states. In particular, gapped QSL 
states exhibit a particularly simple type of long-range entan- 
glement, characterized by Topological Entanglement Entropy 
(TEE) 2627 . Often the two types of characterization are con- 
flated, but this is not necessarily the case. States with both 
long range entanglement, e.g. with TEE, and VBS order exist. 
Such states, while not technically QSLs by the standard defi- 
nition given above, have all the same exotic physics as QSLs 
with unbroken spatial symmetry. We note, however, that it is 
believed that for S=l/2 spins on a lattice such as this one with 
an odd number of spins per unit cell, the absence of VBS order 
implies the presence of long-range entanglement. Therefore a 
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FIG. 1: (Color online) The ground state phase diagram for the spin- 
| AFM Heisenberg J\-J% model on the square lattice, as deter- 
mined by accurate DMRG calculations on long cylinders with L y 
up to 14. Changing the coupling parameter J2/J1, three different 
phases are found: Neel antiferromagnet (AFM), topological quan- 
tum spin liquid (QSL), and stripe AFM phase. m s (ko = (n,n)) 
[m s (k x = (n, 0))] denotes the staggered magnetization in the Neel 
AFM phase [stripe AFM phase], whose saturation value is 1/2. As 
and At denote the spin singlet gap and spin triplet gap, respectively. 



convincing demonstration of vanishing VBS order does, indi- 
rectly, imply interesting QSL physics. It is, however, less im- 
portant to characterizing and proving the existence of a QSL 
than positive, direct evidence of long-range entanglement. 

Most of the literature on the intermediate phase of the J\- 
J 2 model has focused on the possibility of symmetry breaking 
VBS order. Many of these prior studied have suggested that 
the intermediate state has VBS order. We note, however, that 
all numerical results for the J\- J 2 model are based either on 
biased techniques (such as series expansion or coupled clus- 
ter methods, or fixed node or related versions of Monte Carlo 
adapted to avoid the sign problem which is present for unbi- 
ased Monte Carlo in this system), or on exact diagonalization 
of very small systems. Some theoretical motivation for the 
possibility of VBS order comes from the theory of deconfined 
quantum criticality 28 , which predicts that a continuous quan- 
tum phase transition - a deconfined quantum critical point 
(DQCP) - should occur between an ordered Neel state and 
a plaquette or columnar VBS state, in some models. How- 
ever, the existence of such a transition does not in any way 
imply that it occurs for the J1-J2 model in question, or that 
this particular model even harbors a VBS phase. Other theo- 
retical motiviation for VBS order comes from its presence in 
some large-iV generalizations of the nearest-neighbor Heisen- 
berg antiferromagnet. However, these large N studies are not 
controllably close to the SU(2) case and moreover do not con- 
sider second neighbor interactions. In short, we believe there 
is very little compelling evidence for the existence of VBS 
order in the isotropic 5=1/2 J1-J2 model to be found in 
the prior literature. We will return to discuss VBS states in 
Sec. rvTAl 

The only unbiased technique capable of treating generic 
frustrated two dimensional spin systems of moderately large 
size is the Density Matrix Renormalization Group (DMRG) 



method PESO While the sizes that can be studied using the 
DMRG are not as large as those accessibly by quantum Monte 
Carlo (QMC) for unfrustrated models, they are still very large 
and they are not limited by the sign problem, which prevents 
application of QMC to most realistic physical models. More- 
over, the DMRG has some advantages over QMC: it is intrin- 
sically a zero temperature technique, and obtains a convenient 
representation of the ground state wavefunction. Most impor- 
tantly for our purposes, the DMRG is very efficient and conve- 
nient for calculating the entanglement entropy, which we re- 
turn to in some detail below. In this paper, we report the results 
of extensive simulations (with truncation error ~ 1CP 7 ) on nu- 
merous cylinders of circumference L y = 3 — 14, and lengths 
L x > 2L y . In our simulations, we measure spin-spin corre- 
lation functions, correlation functions and expectation values 
of VBS order parameters, bulk singlet and triplet energy gaps, 
and entanglement entropy. All results confirm the existence 
of magnetic order for small and large J2, and that (see Fig.[T| 
the ground state for 0.41 < J2/J1 < 0.62 is non-magnetic, in 
very good agreement with the most accurate prior results from 
series expansion and coupled cluster 24 methods. Furthermore, 
we find that the intermediate phase has a gap to both singlet 
and triplet excitations and, within our uncertainty, no VBS or- 
der in the 2D limit as extrapolated from the VBS correlation 
functions. We carry out further checks for possible finite-size 
effects due to the boundaries, to see if this might artificially 
suppress VBS order, and see no indication that this is the case. 

The latter results suggests a QSL state, based on negative 
evidence: the apparent absence of VBS order. We find two 
positive evidences that this suggestion is correct, and that the 
state is a Z2 QSL. First, we find a non-zero TEE, 7, which 
is a constant and universal reduction of the von Neumann en- 
tanglement entropy, known to vanish in any gapped state with 
short-range entanglement. Notably, we point out in Sec. IV 
that discrete spontaneous symmetry breaking phases such as 
valence bond solids have absolute ground states which are 
Schrodinger cat states with a constant enhancement of the en- 
tanglement entropy - i.e. an effect of opposite sign to the 
TEE. Phases with non-zero 7 and a gap to all excitations are 
topological phases. Like conformal field theories in two di- 
mensions, only discrete types of topological phases exist, with 
discrete allowed values of 7 (which plays a role somewhat 
similar to the central charge in a conformal field theory). For 
all points we have studied within the non-magnetic phase, the 
value of 7 is equal, within numerical uncertainty of 2%, to 
ln(2), which is the minimal value possible for 7 in a topolog- 
ical phase with time-reversal symmetry. A topological entan- 
glement entropy of 7 = ln(2) implies either a Z2 QSL or a 
"doubled semion" phased As there is, to our knowledge, no 
theory suggesting the appearance of the semion phase in an 
SU(2) invariant spin-1/2 model, we take this as strong evi- 
dence for a Z2 QSL state. The second positive evidence for a 
Z 2 QSL is a remarkable odd/even effect, in which static VBS 
order is entirely absent for even L y but is observed directly in 
the VBS expectation values for odd L y . This is expected on 
general theoretical grounds for a Z2 QSL, as we show in Ap- 



pendix A 1 We compare the behavior of the numerically ob- 
served static VBS order for odd circumference cylinders with 
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theory, and find quite consistent results. 

The remainder of the paper is organized as follows. In 
Sec. [n] we report results of magnetic and dimer correlation 
functions, and their extrapolation to the infinite system limit. 



we define the dimer-dimer correlation 
BfB 1 ^), with the corresponding structure factor 



Sec. Ill discusses the singlet and triplet energy gaps. Sec. IV 
describes the theory and measurements of the topological en- 
tanglement entropy, and Sec. [V] presents results on the even- 
odd effect. We conclude in Sec. [Vljwith a summary of the 
conclusions, and a detailed discussion of the reasons to think 
VBS order, even weak, is unlikely in this model, in response 
to a recent critique The Appendix gives a theoretical deriva- 
tion and discussion of some properties of Z2 quantum spin 
liquids. 



n. CORRELATION FUNCTIONS 

In this section we discuss the behavior of correlation func- 
tions of spin and dimer (VBS) operators. Here and in the rest 
of the paper, all our numerical data is based on DMRG simula- 
tions on cylinders, i.e. finite square lattices with N = L x x L y 
sites and with open and periodic boundary conditions in the x 
and y directions, respectively. When not otherwise specified, 
we fix the aspect ratio to L x /L y = 2, with L y = L, then 
L x = 2L, which has been shown to optimize results in the 
DMRCPSED Moreover, to extract bulk properties, we will 
often work on the central half of the system with an effective 
system size N c — L x L. For instance, in computing spin 
correlation functions (Sj • Sj), we restrict site indices i and 
j to the central half of the system so that the obtained corre- 
lation functions could represent the bulk properties. We keep 
more than m = 12000 states in each DMRG block for most 
systems, which is found to give excellent convergence with 
truncation errors of the order or less than 10~ 7 . 

We begin with measurements of the magnetic correlations 
in the ground state, (Si • Sj), and the corresponding static 
structure factor M B (k,L) = ± J2ij e ik (ri ~ rj) (S l • Sj}. The 
structure factor is peaked at k = (tt, tt) for small J 2 and 
k x = (tt,0) or k y = (0, tt) for large J2, corresponding to 
the Neel and striped AFM states, respectively. To quantita- 
tively analyze the order, we perform an extrapolation of the 
(squared) staggered magnetization, m 2 s {k, L) = -JaM s (k, L), 
to the two dimensional limit (L = 00) according to the gen- 
erally accepted form m 2 s (k,L) — m? s (k, 00) + j- + (see 
Figsfga) and (b)). 

Extrapolation from data for L < 12 shows that the Neel 
AFM order is non-zero for J2 < 0.41, while striped AFM or- 
der onsets for J 2 > 0.62, thus establishing the phase bound- 
aries shown in Fig.[T] A strong check on the quality of our re- 
sults is the staggered magnetization at J2 = 0, which we find 
to be m s (ko, oo)=0.304, very close to the best known numer- 
ical value of the magnetic moment m s =0.307 by large-scale 
quantum Monte -Carlo (QMC) simulation 33 . The locat ion o f 
the phase boundaries is consistent with previous studies^H. 

We next consider possible VBS order, which has been 
considered a prime candidate for non-magnetic symme- 
try breaking in the intermediate phase. From the bond 
operators Bf = S; • Si+ Q on bond (i, i + a) with 



a — x 
functions 

Mf(k,L) = j2j2ij elk '^ l ~ rj ^ {{B?Bj) - 
Typical VBS patterns expected theoretically have momentum 
k x = (ir,0) or k y — (0,tt), so to study the correlations, we 
focus on L y even, for which k y = tt is an allowed momen- 
tum. We indeed observe a maximum in M^ a (k, L) at k = k a 
(a = x, y), and therefore define the dimer order parameters by 



j2 M% a (k a . L). As shown in the inset of Fig 



for finite systems, both horizontal and vertical dimer ore 
parameters have a maximum within the intermediate phase. 
Note that for the larger systems, the order parameters for hor- 
izontal and vertical dimers become nearly indistinguishable, 
indicating that the isotropy of the two dimensional limit is be- 
ing recovered. 

Applying the same extrapolation scheme used for the mag- 
netic order parameters, however, the extrapolated dimeriza- 
tion m 2 , a (see Figj4 a)) for L — >• 00 vanishes for all < 
J2 < 1. For characteristic values of J2 near the middle of 
the intermediate phase, an exponential fit of the dimer-dimer 
correlation function (not shown) gives an estimate of the VBS 
correlation length £4 rj 4. Taken at face value, these obser- 
vations indicate that the VBS order is a finite-size effect, and 
vanishes in the thermodynamic limit. More conservatively, at 
a minimum, the result indicates that the VBS correlations we 
observe cannot be distinguished from just fluctuation effects 
in a state with unbroken spatial symmetry, and there is no a 
priori reason to regard them as evidence of true VBS order. 

Both columnar and plaquette VBS phases have been sug- 
gested in the past. The complex order parameter rn^.x +imd,y 
in fact is sufficient to detect and distinguish both columnar 
and plaquette VBS phases 28 , but as an additional check we 
measure directly the correlations of the plaquette operator 
Pi = ^(Ili + II^ 1 ) where II; cyclically permutes the four 
spins of the plaquette i in a clockwise fashion. The plaquette 
order parameter determined from the corresponding structure 
factor (see Supplementary Information) is shown in Fig. |4|b). 
Like the VBS order parameter, it vanishes in the extrapolation 
to the thermodynamic limit. 



III. ENERGY GAPS 

We next consider the energy gap to bulk singlet and triplet 
excited states, and find both to be non-zero in the intermediate 
phase. This rules out any type of magnetic order, not just the 
(tt, tt) and (tt, 0) orders considered explicitly via the correla- 
tion functions. It also rules out other exotic states breaking 
SU(2) symmetry, such as spin nematics. This is because any 
state with broken spin-rotational symmetry must have a van- 
ishing gap by Goldstone's theorem. 

To obtain bulk excited states, we follow Refs. 7 30 31] and 
first target only one state, sweeping enough to obtain a high- 
accuracy ground state; then we restrict the range of bonds 
that are updated in the DMRG sweeps to the central half 
of the sample and target the two lowest-energy states, again 
sweeping to high accuracy, but keeping the end regions of 
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FIG. 2: (Color online) Finite-size extrapolations of the magnetic or- 
der parameters and spin excitation gaps, (a) The Neel AFM order 
parameter m;?(k) at wavevector ko = (tt, 7t) and (b) stripe AFM or- 
der parameter m^(k) at wavevector k^ = (tt, 0) or k^ = (0, tt), for 
various values of J2, fitted using second-order polynomials in 1/L. 
Neel AFM order disappears for J2 > 0.41, while stripe AFM or- 
der develops for J2 > 0.62, as seen in the corresponding insets, (c) 
Spin triplet gap At and (d) spin singlet gap As for different values 
of J2, also fitted using second-order polynomials in 1/L. The inset 
in (c) shows At for L = 4, 6, 8, 10, and the extrapolated values in 
the 2D limit, as functions of J2 . For the spin singlet gap, due to the 
numerical cost, we focus on several typical data points as shown in 
(d). 
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FIG. 3: (Color online) Order parameters for horizontal and vertical 
dimers. (a) The dimer order parameter rr? d x at wavevector k x = 
(tt, 0) and (b) plaquette order parameter m^ y at k v = (0, tt), as a 
function of J 2 for different system sizes, and extrapolated to L = 00. 
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FIG. 4: (Color online) Finite-size extrapolations of the dimer order 
parameter and plaquette order parameter, (a) The dimer order pa- 
rameter rri^ y at wavevector k H = (0, tt) and (b) plaquette order 
parameter trip, for various values of J2, fitted using second-order 
polynomials in 1 / L. The inset shows the plaquette order parameter 
for L — 4,6, 8, 10, and the extrapolated values in the 2D limit, as 
functions of J2 . 



the samples locally in the ground state. To obtain the spin 
triplet gap, we do similar things, but target states with total 
S z = and S z = 1 separately. As for the staggered mag- 
netization, we perform a second order polynomial extrapola- 
tion of the singlet and triplet gaps to the thermodynamic limit 
(Figs.|2|c,d)). Consistent with expectation, both A$(L = 00) 
and At(L — 00) vanish in the two AFM phases. They are 



5 



both, however, non-zero and large in the intervening region 
(see Fig. [T}. This rules out any state with broken SU(2) spin 
symmetry. 

We notice that the singlet gap remains consistently below 
the triplet gap throughout the intermediate phase. This is an 
indication of short-range singlet formation. It is consistent 
with a spin liquid state, and with a system with weak VBS or- 
der. We would, however, expect a strong VBS state to have a 
triplon excitation, corresponding to breaking one singlet bond, 
as the lowest energy bulk excitation, lower than singlet excita- 
tions which require breaking two singlets. So we can exclude 
a strong VBS state in this sense based on the excitation spec- 
trum. 



IV. TOPOLOGICAL ENTANGLEMENT ENTROPY 

The above results provide evidence against conventional or- 
dering in the intermediate region. Magnetic ordering appears 
comfortably excluded by both the correlation function and ex- 
citation spectrum analysis. Extrapolation of the dimer corre- 
lations to the thermodynamic limit argues that VBS order is 
absent as well, but we cannot exclude some very weak order- 
ing on these grounds alone. 

We now undertake a positive evidence for a QSL with topo- 
logical order - the topological entanglement entropy. The 
topological entanglement entropy is obtained from the von 
Neumann entanglement entropy S(A). The latter is defined 
for a state |^>o) (which we take to be the ground state) and 
a partition of the full system into a subsystem A and its 
complement B, by first constructing the reduced density ma- 
trix pa = Trs| , 0o)('0o|- Then the entanglement entropy 
S(A) = — Tr a(pa In pa)- For a system with a gap to all 



bulk excitations (as we have verified in Sec. Ill i, provided the 
boundary between A and B is taken to be smooth (i.e. have 
no corners), the entanglement entropy must scale according to 



S{A) ~(tL-j- 



(2) 



where the omitted terms vanish in the large L limit. Here a is 
a non-universal number that measures the local entanglement 
across the boundary. According to Refs. 1261271 the po sitive 
term 7 is universal constant reduction from the area lawP^l 
It arises entirely from non-local entanglement, and is topolog- 
ical in origin. In particular, the area law is strictly obeyed, i.e. 
7 = 0, for any state without long-range entanglement, that is, 
which can be smoothly deformed into a product state. This 
is true, in the absence of spontaneously broken symmetry, for 
any ground state which does not exhibit topological order, i.e. 
which is not a topological QSLP^ZI 

Although it is not discussed in the seminal papers on topo- 
logical entanglement entropy, a non-zero negative 7 (i.e. a 
positive correction to the area law) can arise from discrete 
spontaneous symmetry breaking (more severe positive correc- 
tions to the area law arise in the case of a continuous broken 
symmetry) 35 ! 36 ! b u t this is inconsistent with the existence of a 
gap to bulk excitations). In particular, in an ideal model with 
an exact discrete symmetry of the Hamiltonian, the eigen- 
states must form irreducible representations of the symmetry 



group. For simple abelian groups such as Zn, these represen- 
tations are one dimensional, so this implies the Hamiltonian 
eigenstates are mutual eigenstates of the symmetry genera- 
tors. This applies of course to the absolute ground state of 
the system, which is therefore a Schrodinger cat state, which 
superimposes the symmetry broken global ground states with 
equal weight. For the case of a fully broken Z^ symmetry, 
with N degenerate ground states in the thermodynamic limit, 
this gives rise to N terms in the Schmidt decomposition of 
the ground state, and therefore a correction 7 = — In (AT), 
i.e. a positive correction to the area law or ln(N). We have 
indeed observed such behavior numerically in test studies of 
the simplest quantum transverse field Ising model in the fer- 
romagnetic phase, consistent with the expected 7 = — ln(2) 
for this case. 

Thus we see that there are two potential sources of a non- 
zero constant term in the entanglement entropy. A topologi- 
cal contribution which decreases the entropy, and a symmetry 
breaking contribution which increases it. The latter correction 
arises from global entanglement of the entire system. In work 
completed since the earlier version of this article appeared,^ 
it has been shown that the DMRG, which is a minimum entan- 
glement approximation, tends to converge, for large systems, 
to quasi-ground states which capture all entanglement out to 
a long length scale, but not the last global entanglement. That 
is, for long systems, the convergence of the DMRG is first to 
a Minimum Entanglement State (MES) amongst the manifold 
of states comprising the degenerate ground states in the ther- 
modynamic limit 37 . For topologically ordered phases, which 
have a ground state degeneracy in the thermodynamic limit of 
topological origin, the MES exhibits the universal reduction 
of entanglement entropy, i.e. the universal positive value of 7. 
For symmetry broken states, for which there is a ground state 
degeneracy in the thermodynamic limit dictated by symmetry, 
the MES is simply a single product-like state, with 7 = 0. As 
shown by Jiang, Wang and Balents in Ref.[37] for a fixed sys- 
tem size which is not too large, the DMRG can be pushed to 
converge to the global ground state, by increasing the number 
of states m. This is accompanied in these cases by a sharp in- 
crease in the entanglement entropy. By increasing the length 
L x of the system at fixed L y , this final increase in the en- 
tanglement entropy can be pushed beyond the range of fea- 
sible calculations, and the simulation is guaranteed to obtain 
the MES. In the MES, the constant correction 7 is entirely of 
topological origin, and is zero in discrete symmetry breaking 
states. Thus in this limit 7 is the topological entanglement en- 
tropy, and a non-zero result proves that the state is a (topolog- 
ical) QSL. Moreover, we see from the above discussion that 
a positive 7 can only come from topological order, so we do 
not obtain false positive signatures of topological order from 
symmetry breaking. 

In Figj5ja), we plot von Neumann entanglement entropy 
S(L y ) associated with the constant x cut which separates the 
cylinder into two symmetric parts of equal length, L x /2, as a 
function of L y , with L y even (for L y odd, there are additional 
effects which we discuss in Sec. [Vj. By comparing systems 
of different lengths (Fig.[3J>), we see that the entropy is essen- 
tially independent of L x for L x > 2L y , and so equal to its 
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with (even) L x — > oo and odd L y , the Z2 QSL induces a non- 
vanishing staggered dimerization, 
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FIG. 5: (Color online) The entanglement entropy at J2 = 0.5 and 
0.56. (a) The entanglement entropy S(L y ) for L y = 6 — 14. By 
fitting S(Ly) = aL y - 7, we obtain 7 ~ 0.70 ± 0.02 for J 2 = 0.5, 
and 7 ~ 0.72 ± 0.04 for J 2 = 0.56. (b) Length dependence of 
the entanglement entropy for J2 = 0.50 and several system widths. 
One observes that the entropy is almost independent of L x for long 
systems (a small increase with L x can be observed for the smallest 
L x at L v = 8). 



limit at L x = 00. We then extrapolate 7 from the fitting func- 
tion S(L y ) = aLy—j. For J2 = 0.5, deep in the magnetically 
disordered phase, our results show that 7 = 0.70 ± 0.02. This 
value appears constant, within numerical uncertainty, within 
the intermediate phase: for J 2 — 0.56 (close to the quantum 
phase transition point J 2 = 0.62), we obtain in the same way 
7 = 0.72 ± 0.04. Without even consider the magnitude of 7, 
the fact that we see a negative rather than positive correction 
to the entanglement entropy is strong evidence against VBS 
order. 

As mentioned in the Introduction, the topological entangle- 
ment entropy 7 takes discrete values in topological phases. 
The minimum possible value for systems with unbroken time- 
reversal symmetry is 7 = ln(2) ps 0.69, which is within 
2% of the numerical results. The constancy of the numeri- 
cal topological entanglement entropy and the consistency with 
the theoretically allowed value of ln(2) constitute strong evi- 
dence for a topological QSL state. The appearance of the pure 
number ln(2) (within happily small numerical uncertainty, of 
course) is certainly very striking coming out entirely unso- 
licited from the DMRG calculations. 

Notably, 7 = ln(2) is the expected value for a Z2 QSL 
phase. The Z2 QSL is in many ways the simplest spin-liquid 
state, and has appeared repeatedly in theories of quantum 
magnets. As a rather complete theory of the low energy prop- 
erties of Z2 QSLs is available, we can compare this to numer- 
ics in various ways. 



V. ODD-EVEN EFFECT 



(£?>=B* + l> B (-i)*s 



(3) 



with D, 



-Ly/S 



exponentially decreasing with circumfer- 



ence. By contrast, no dimerization appears for even L y . We 



obtain this behavior in Appendix A 1 directly from the effec 



tive Z2 gauge theory description, which shows that it is a uni- 
versal feature of Z2 QSLs on the square lattice, and not par- 
ticular to the quantum dimer models studied in Refill 

Precisely this behavior is observed in our numerics. 
Fig. [6ja,b) contrast the oscillatory and non-oscillatory hori- 
zontal bond expectation values obtained for odd and even L y . 
For even L y , some small boundary effects are observed, de- 
caying over ~ 3 lattice spacings. Figj6jc) shows the exponen- 
tial behavior of D x obtained as the difference of even and odd 
bonds at the center of the sample. Interestingly, theories pre- 
dict (see Ref 38 and also Supplementary Information) £ = 2£, 
where £ is the true dimer correlation length defined through 
the dimer correlation function. This explains the rather slow 
decay of D x , which fits to £ ps 5, reasonably consistent with 
£d fa 4 found (see Sec. |ll| from the examination of VBS cor- 
relation functions. While some even-odd effect might be ex- 
pected in a columnar dimer phase for narrow cylinders, the 
exponentially-decaying behavior and results of other tests (see 
Sec. |VI A| > seem consistent only with a Z2 QSL. 



VI. DISCUSSION 

The previous sections have shown that DMRG makes a 
compelling case for a non-magnetic intermediate state in the 
J\-J 2 model. From direct measurements of the dimer order 
parameter and correlations, the intermediate state appears to 
have no or very weak VBS order. Most dramatically, we find 
a robust constant suppression of the entanglement entropy rel- 
ative to the generic area law, known as topological entangle- 
ment entropy, which is a unequivocal signature of topological 
order. The value of the topological entanglement entropy we 
find is within 2% (and our numerical uncertainty) of the ex- 
pected universal value 7 = ln(2) for the simplest Z 2 QSL 
state, which suggests comparison of specific theoretical pre- 
diction for this Z2 phase to numerics. We indeed find a char- 
acteristic even-odd effect in the staggered dimerization, con- 
sistent with this state. 

It is worth noting that ours is not the only suggestion of a 
QSL state in the J\-J% model. Notably, after the initial ver- 
sion of this paper appeared, a parallel wortP^came to similar 
conclusions based on a tensor network variational method. 



In this section, we make such a comparison based on the 
theory of the Z 2 QSL. Specifically, in a Z 2 QSL on the square 
lattice, it is predicted that cylinders with odd circumference - 
and not those with even circumference - should exhibit non- 
vanishing bulk staggered dimerization. This even-odd effect 
was first obtained, to our knowledge in Refl38l by analysis of 
quantum dimer models^2H2l Specifically, for a long cylinder 



A. Could this be a weak VBS state with strong finite size 
effects? 

In our opinion the above results all point in the same direc- 
tion, and are especially definitive given the seemingly unas- 
sailable implication of the observed topological entanglement 
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entropy. Nevertheless, following an earlier version of this pa- 
per, SandvilipJ has suggested, by comparison with quantum 
Monte Carlo results for so-called J-Q models on cylinders, 
that similar behavior might occur for a system with a VBS 
ground state in the thermodynamic limit, due to strong finite 
size effects. We discuss this suggestion here. 

1. Difference of models 

The results of Ref.l32lare based on the J-Q models, which 
have four or six spin interactions (with coefficient Q). These 
multi-spin interactions explicitly involve interactions between 
dimers, and as a consequence rather naturally favor VBS 
states. For instance, the simplest mean-field treatment of the Q 
term in the J-Q2 model would proceed from by decoupling it 
by defining a mean-field dimer expectation value of the dimer 
operator, and thereby a VBS phase appears when the Q term 
becomes substantial. Thus it is natural and intuitive to expect 
a VBS phase in the J-Q models. By contrast, there is no a 
priori reason to expect dimer order in the J1-J2 model. The 
notion that a VBS state is somehow the most "likely" candi- 
date for the intermediate non-magnetic state in the 3\-J-x case 
is a misleading starting point. More importantly, we should be 
cautious in drawing conclusions from the J-Q models on the 
behavior of the J\-Ji model. 



2. Entanglement entropy 

The most direct evidence for a QSL state we have obtained 
is the topological entanglement entropy, remarkably close to 
the universal expected value for a Z2 QSL. In Ref. 32, Sand- 
vik suggests that "it would not be surprising" if a system near 
a Neel to VBS transition (i.e. a DQCP) would exhibit a con- 
stant correction to the area law similar to that expected for a 
topological phase. In fact, we have shown theoretically that 
in a VBS state, there is indeed a constant correction but of 
opposite sign to that of a topological phase. Thus even forget- 
ting the magnitude of 7, the sign alone is a strong argument 
against VBS order. The fact that the measured 7 is within 
2% of the very beautiful and universal expected result ln(2) 
makes it hard to imagine this is mere coincidence. 

In the context of a putative DQCP, the constant correction 
for a VBS state obtains if the system size is larger than the "de- 
confinement length", below which there is an emergent U(l) 
symmetry unifying the plaquette and columnar VBS states, 
and linear combinations in between, with one another. Would 
one perhaps see a signal similar to the topological entangle- 
ment entropy were this length longer than the system size? 
Actually in this case we expect that the system should ap- 
pear to exhibit a gapless Goldstone mode, characteristic of 
spontaneously breaking this U(l) symmetry. This is the sit- 
uation discussed in Refs.35 36 In fact the behavior of the 
entanglement entropy in this case is even further from that of 
a topological phase: a positive logarithmic enhancement of 
the entanglement entropy beyond the area law is predicted, 
again of opposite sign to the topological case. Moreover, in 



the Id limit, L x ^> L y , the system should behave as a 1+1- 
dimensional conformal field theory with central charge c = 1, 
and hence exhibit a logarithmic growth of entanglement en- 
tropy, S(L) ~ | ln(L x ), in such a case. This is completely 
at odds with our observations - observe the constant behavior 
versus L x in Fig. [5j?. 

3. VBS scaling 

In Fig. 23 of Ref 32 our data for the dimerization is replot- 
ted along with data for the J-Q2 model on a log-log plot, to 
fit to a single pure power law. Data for g(= J2/J1) = 0.5 is 
compared to Dy ~ L~ a with a m 1.8, and for g — 0.56 is 
slightly above it. Small details of the data for the latter case for 
the smallest systems, L x = 4, 6, are used to conclude that the 
system is VBS ordered in the infinite-size limit. We disagree. 
First, note the simple fact that the dimerization for the J1-J2 
model is much smaller than that of the J-Q2 model (which is 
the more weakly VBS ordered of the two J-Q models). Sec- 
ond, the scaling on this plot for the J1-J2 model (unlike the 
J-Q models) is quite close to a = 2, which is, as mentioned in 
the same paragraph of Ref. 32 exactly the behavior expected 
for a non- VBS phase. 



4. Even-odd effects in VBS states 

One of the pieces of evidence for the Z2 QSL state taken 
from our numerics was the very distinct behavior of the stag- 
gered dimerization in even and odd circumference systems, 
described in Sec. [V] While this is certainly consistent with a 
Z2 state, one could imagine similar behavior arising in a sys- 
tem with VBS order in the thermodynamic limit. Here we 
consider the expected behavior in such a situation more care- 
fully, for comparison to our results. 

Consider a system which is spontaneously dimerized in the 
2d limit, with a columnar dimer ground state. This state is 
four-fold degenerate, with four ground states consisting of two 
states with "horizontal dimers" staggered along the x direction 
((■7T, 0) order), and two states with "vertical dimers" staggered 
along the y direction ((0, ir) order). In the thermodynamic 
limit, these states are degenerate by rotation and translation 
symmetry. When confined to a cylinder, the anisotropy of the 
boundary conditions breaks the symmetry between the hori- 
zontal and vertical states. For the case of odd- width cylin- 
ders, the vertical dimerization is frustrated, because alternat- 
ing "rows" of vertical dimers do not fit into the sample. This 
clearly would favor the horizontal dimer states. Amongst the 
two horizontal dimer states, the presence of an end to the sys- 
tem splits the remaining degeneracy, so all degeneracy is bro- 
ken and we would expect long-range horizontal dimer order 
to appear. To this extent, the behavior for odd-width cylinders 
is the same as observed in our numerics, and as expected for 
the Z2 QSL. The difference is in the scaling. If the 2d system 
has a gapped dimer ground state, we would expect the expec- 
tation value of the dimerization to converge exponentially to 
a non-zero two dimensional limit as the width of the cylinder 
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FIG. 6: (color online) Even-odd effect, (a) Expectation value of 
horizontal bond operator, (Bf), for L y = 5, L x = 32. (b) The 
same expectation value for L y — 6, L x = 32. (c) Dimer order 
parameter Dd, x for odd L y at L x = oo. The red line denotes the 
exponential-decaying fitting function with the form in Eq.j4]l. (d) 
Modified boundary induced dimer order parameter for Ly = 6, 8, 
with d the distance from the boundary. Here the dimer order param- 
eter is defined as the dimer density difference between two nearest 
neighbor vertical dimer bonds. Inset shows the correlation length £ 
along the cylinder as a function of L y . 



increases, i.e. 



D 



x 1 2d dimer state 



(4) 



where A and £ are constants, and Drx, is the value of the dimer 
order parameter in the thermodynamic limit. 

As shown in Fig(6jc), the numerical fitting to this form 
gives D oo = within numerical accuracy. This is en- 
tirely consistent with vanishing VBS order and a Z^ QSL in 
the thermodynamic limit, but of course cannot exclude some 
weak dimerization smaller than our numerical uncertainties. 
It seems to us natural to take the former interpretation, since it 
is simpler. Accor ding to the theory for the Z2 QSL discussed 
in Appendix A 1 and quantum dimer model results 38 , an expo- 
nential decay for the staggered dimerization is expected with 
"doubled" correlation length £ = 2£, where £ is the true VBS 
correlation length. We find £ ss 10 lattice spacings, which is 
equivalent via Eq. ( |A28| > to £ « 5. 

For the even circumference cylinders, the vertical dimer or- 
der is unfrustrated, and it is an energetic question, which likely 
depends upon the details of the model, whether the vertical or 
horizontal dimer order would be favored in this case. If the 
horizontal dimer state is favored, then we again expect be- 
havior like Eq. Q, which is manifestly inconsistent with our 
numerics, and markedly different from the Z2 QSL. However, 
it is perfectly conceivable that the vertical dimer pattern is fa- 
vored instead. If so, the periodic boundary conditions do not 
break the symmetry between the two vertical dimer states, and 
so we expect the DMRG to converge to the symmetric linear 
combination of the two dimer states, which lacks any sponta- 
neous dimer pattern. So at least the presence of an even-odd 



effect in the static dimerization is consistent with a VBS state, 
if the cylindrical geometry favors the two VBS states with hor- 
izontal rows of vertical dimers. On the face of it, this appears 
consistent with our numerical results for the staggered dimer- 
ization, if one assumes that the value of the dimerization itself 
(extrapolated from odd circumference cylinders) is smaller 
than our numerical uncertainty. But it is worth pointing out 
that for this scenario to hold, the even circumference system 
must be in a Schrodinger cat state, and should exhibit a posi- 
tive ln(2) enhancement of the entanglement entropy (negative 
TEE) as a consequence, and moreover convergence to such 
a state should be progressively more difficult with increasing 
L T . This is not at all what we see. 



5. End effects 

In Ref. |32l strong boundary effects are observed on the 
dimerization in the J-Q models. Indeed, on symmetry 
grounds, an open end breaks translation and reflection sym- 
metries in the x direction, and as such should act as a "bound- 
ary field" on the staggered dimer order D x , i.e. it induces a 
term — \D x (x = 0) in a Landau theory of this order. On these 
grounds, we always expect some staggered dimer order near 
the boundary. If it is energetically disfavored in the bulk, this 
will decay rapidly. Otherwise, it will penetrate deep into the 
bulk. In the J-Q models, it was found that the boundaries in- 
duce a quite strong dimerization, so that for even L y the bond 
expectation values (Bf) oscillate visibly (c.f. in the inset of 
Fig. 6, and in Fig. 15a of Refl32l the bond expectation value 
shows oscillations with large amplitude in the J-Q 3 and J-Q 2 
models, respectively). By contrast, in the J1-J2 model, we 
see in Fig. |6|b) that there are no visible oscillations in the 
same quantity when L y is even. This qualitative difference 
tells us that D x order is clearly much less favorable in the Ji- 
J2 model. 

We next try to address the possibility, raised above, that the 
cylindrical geometry when L y is even favors D y order, i.e. 
horizontal rows of vertical dimers. This is at odds with our 
measurements of the dimer correlations and the entanglement 
entropy. Still, it is more compelling to explicitly test to rule 
out the possibility directly. To do so, we have studied sev- 
eral modified cylinders with even circumference, in which the 
ends of the cylinder have been altered, breaking translational 
symmetry along y in order to break the degeneracy and fa- 
vor one of the two vertical dimer states. What we observe is 
that in all cases, as shown in Figj6fd), although dimer order 
is induced by this symmetry breaking in the vicinity of the 
boundary, it decays exponentially into the bulk of the cylin- 
der. The correlation length £„ for this vertical dimer order still 
depends on circumference for the system sizes in our study, so 
we plot it versus L y to see if it is limited by the system size (it 
does not appear to be), and to extrapolate from this its value 
in the thermodynamic limit. We observe that this correlation 
length grows sub-linearly in L y , and extrapolates to £„ ~ 4 
in the 2D limit (i.e., L y = 00). This is very different from 
what would be expected for a 2d state with long-range dimer 
order, in which the non-zero stiffness (surface tension) of the 
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ordered dimer state would prevent such decay (we would ex- 
pect = oo in this case). If one were to imagine that the sys- 
tem were proximate to a DQCP, and L y were smaller than the 
deconfinement length, then we would instead expect £ v oc L y , 
which again is not consistent with our results. Note also that 
the value for is quite consistent with the value for £ ob- 
tained earlier. The fact that vertical dimer order decays away, 
even when the most favorable conditions have been created 
for it, is strong evidence against VBS order in the 2d limit. 

B. Summary and Open Issues 

In conclusion, we have presented compelling evidence from 
accurate DMRG calculations for a topological QSL state in 
the two dimensional J\-Ji Heisenberg model. This is the sim- 
plest example of such a QSL discovered to date, and the only 
one to our knowledge for a Heisenberg model on a Bravais 
lattice. As such, it is particularly attractive for further theoret- 
ical and experimental study. We anticipate, for instance, that 
our discovery will afford an opportunity to explore the QSL 
mechanism of unconventional superconductivity^^ in a con- 
trolled theoretical setting. 

Another consequence of topological order is the presence 
of quasi-degenerate ground states on the torus or cylinder. A 
two-fold quasi-degeneracy is expected for a Z-i QSL on the 
cylinder studied here, with a splitting of order L x e~ Ly ^ in 
the case of long cylinders, where £ is the spin-spin correlation 
length (see Appendix |A2| i. As shown in Refl37land discussed 
in Sec. |IV| the DMRG preferentially converges, however, to 
just one of the quasi-degenerate ground states (specifically, 
a minimally entangled state). This explains the absence of 
an observed topological degeneracy in this and other DMRG 
studiesPESSI it is a non-trivial and open problem to obtain 
the second ground state and thereby extract the topological 
energy splitting. It is our expectation that it is actually orders 



of magnitude smaller than the bulk energy gaps. 

The nature of the quantum phase transitions from the QSL 
to Neel and striped antiferromagnetic phases is an interesting 
topic for future study. Though we have not focused on the 
transitions themselves, and more work is clearly required to 
make strong conclusions about them numerically, it appears 
that the transition from the Neel to QSL state may be continu- 
ous. Ref (32] erroneously claims that a Neel to QSL transition 
might be in the same universality class as the DQCP between 
Neel and VBS order, because "the operator causing the VBS 
order is dangerously invariant". Though at the DQCP the op- 
erator which distinguishes between columnar and plaquette 
VBS order is dangerously irrelevant, even when this opera- 
tor's coefficient in the Hamiltonian is tuned to zero, the non- 
magnetic phase has spontaneous VBS order. So this claim is 
incorrect. In fact, such a transition requires an entirely differ- 
ent theory. A novel suggestion for the theory of this critical 
point has been made in Ref.44j and it would be interesting to 
compare it to further numerical studies. 
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Appendix A: Z 2 gauge theory 

Here we discuss an effective Z 2 gauge theory description 
of the QSL state, 45 and in particular derive the behavior of the 
dimerization and ground state quasi-degeneracy discussed in 
the main text. We begin with the Hamiltonian 



h = - k h n < >>y,< 

a {ij)ea {ij 



(Al) 



E 

(ij) 



<j[ 



where n, 



= b L b ia 



tb\ a b ja 



and r\i 



(haCapbja + h.C.) 



(-1)* 



We introduced 



"spinon" operators bi a which transform as spinors under 

SU(2), and obey standard commutation relations [b ia , &t„] = 
§ . ..... . ..■<-. ■ 



ij5 a p. The physical spin operators are related to them by 



hbl 



Papb, 



*J8- 



The of operators are Pauli matrix Z 2 



gauge fields, which we will refer to as the "magnetic" gauge 
fields. Z 2 gauge symmetry is enforced by the constraint 



n - 

b'-*l=i 



-(-1)' 



Note that the product in Eq. ( A2 1 is over j not i. 
analog of Gauss' law for the "electric" field afj. 
straint "generates" the Ising gauge symmetry of^ 



(A2) 

This is the 
This con- 

-> s i s j°ij> 



b- t — > Sibi, where Si = ±1 can be chosen arbitrarily for each 
site. 



1. Staggered dimerization 

Here we obtain the behavior of the staggered dimerization 
from the Z 2 gauge theory. For this purpose, it is sufficient to 
integrate out the spinons, since we discuss local properties of 



the QSL state which has a spin gap (but see below Sec. A2i. 
We can obtain this limit from Eq. (All by taking r large, 



which projects the problem onto the subspace with m = 0. 
Then the Hamiltonian reduces to 



H = - K J2 II <i >'Y.<r 



and 



n 

|J'-»I=1 



a ij = 



(A3) 



(A4) 



Eqs. ( A3|A4 i describe the "odd Ising gauge theory". It is in 
the deconfined (QSL) phase for K/h > x c , where x c is some 
order one number specifying the critical point. 

Now consider the staggered dimerization, D x — 
(—l) Xi (Df) — (Df +X ), defined in the main text. On sym- 
metry grounds, we expect that (Df ) oc (af (this relation 
can also be derived by perturbation theory in t/r). We will 
derive the odd/even effect for the staggered dimerization in 
finite-width cylinders in two ways. First, we obtain it directly 
from the Ising gauge theory in the strong coupling limit, which 
is a very short derivation. Second, we obtain it using dual- 
ity and field theory, which exposes the universal nature of the 
staggered dimerization and its relation to Z2 vortex ("vison") 
excitations. 

To see how one might expect the dimerization, we first con- 
sider the "topological" operator 



y=l 



(A5) 



This operator commutes with H and is thus a constant of the 
motion. Moreover, if we consider the case x = 1 at the left 
hand side of the system, we obtain 



Qi 



n n - 

v= l \b'-i|=i 



(-1) J 



(A6) 



where we have used Eq. (A4i. Again using Eq. (A4i, one 
obtains 



(-iy 



(Al) 



Thus Q x = 1 for even L y , but oscillates, Q x = {— \) x , for 
odd L y . Although this is not the dimerization itself, it suggests 
the presence of staggered dimerization in the case of odd L y . 
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a. Direct derivation 

We now turn to the first derivation, working deep in the 
deconfined phase, taking K 3> h, and proceed by direct cal- 
culation perturbatively in h. For h = 0, the ground state(s) are 
obtained by simply choosing a classical configuration of crfj 
with zero Ising gauge flux, fLy^en a ij — 1 on a H plaquettes 
(for instance the state with of.- = +1 on all bonds), and then 
projecting this state to satisfy Eq. \AA\ : 



= i 



where 



P = 



\ n < 

ij-»i=i 



(A8) 



(A9) 



In this state, the expectation value of afj vanishes. This can 
be seen as follows. Define the Wilson loop operator 



w[c] = n < 



(A10) 



where C is a closed curve on the lattice. All such Wilson 
loops commute with the projectors Pj, so \ipo) is an eigen- 
state of the Wilson loop with W[C] IV'o) — IV'o)- Moreover, 
since W[C] 2 = 1, we have 

(VoK#o) = (rk\W[C\ afj W[C}\t/j ) = -(Vok?#o) = 0, 

(All) 

if we choose C to be a curve containing the both (ij). To 
achieve a non-zero result, we must consider non-zero orders 
of perturbation theory in h/K. In general, the form of the 
perturbative eigenstate is 



oo 



\RH' 



IV>o>, 



(A12) 



where R = P(E — H )~ 1 P is the resolvent with H a 
II { h = 0) and Bq the ground state energy of Hq) and 
P = 1 — | ipo) (ipQ | is the projector onto the unperturbed excited 
state subspace, H' = H — Hq = — ^X^(y) a ij> an d trie c„ are 
numerical coefficients. This can be expanded to give a series 
of terms, each involving a product of n electric gauge fields 
acting on |^o) at 0[(h/K) n ]. For each such term, we can re- 
peat the argument in Eq. (jATTi. We will achieve a vanishing 



result provided we can choose C to contain an odd number of 
links that coincide with the set of links C containing the elec- 
tric fields in the corresponding term in the wavefunction and 
the link (ij) in the expectation value. This is always possi- 
ble unless the "dual" of C forms a closed loop. This dual is 
formed by associating a link of the dual lattice with each link 
in C. If the dual of C indeed forms a closed loop, then the 
closed loop C must intersect it an even number of times. 

Thus we obtain non-zero contributions only from terms in 
which C is comprised of closed dual loops. There are triv- 
ial contributions from short loops, the minimal one being the 



case when C contains (ij) twice, which is first order in h/K. 
This gives a non-zero constant contribution to the expectation 
value, but one which is uniform, and hence does not corre- 
spond to a staggered dimerization. A non-trivial result is ob- 
tained first at 0[(h/ K) Ly ^ 1 }, from the smallest closed dual 
loop encircling the cylinder and containing the bond due to 
(ij), which must be a horizontal bond. This leading term 
arises from the 0[(h / K) m ] correction to the ground state ket 
and the 0[(h/ K) Ly ~ 1 ~ m ] correction to the ground state bra 



(m = 0,1,- 



1), giving 
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Ly ~ 1 fl - i\ / h \ t^- 1 ) 



ni—0 



m I \K 



= ■■■ + c n 



L v -1 



( — \) L v x _ 



Here C n is a numerical coefficient which should be deter- 
mined from a more refined analysis. We therefore conclude 
that for odd L y , we obtain the staggered dimerization dis- 
cussed in the main text, with amplitude D x ~ (2h/ K) L «~ 1 = 
exp[—\n(K/2h){L y — 1)], exponentially decaying with cir- 
cumference as advertised. This result derived from the odd 
Ising gauge theory is qualitatively consistent with the one ob- 
tained from the analysis 38 of quantum dimer models. 



b. Dual derivation 

While the above derivation is simple and direct, it relies 
on the strong coupling expansion, which, although it is ex- 
pected to be qualitatively correct in the deconfined phase, is 
not obviously general. It is instructive to obtain the staggered 
dimerization by a more circuitous dual route, which exposes 
the universality of the result and gives a more direct physical 
picture. 

The duality transformation of Eqs. ( A3|A4 i is accom- 
plished by defining 



n 



(ij)ea 



(A14) 
(A15) 



where r Q are new Pauli matrices. In Eq. (A14i (ij) are the 



bonds associated with dual site a at the center of a direct 
plaquette, and in Eq. (A15i, the dual sites a, b are those at 



the centers of the two plaquettes neighboring the bond (ij). 
The scalars u, a b must be chosed to satisfy Eq. ( |A4| i, which re- 
quires that their product around a dual plaquette must equal 
— 1. The dual Hamiltonian is then a fully frustrated transverse 
field Ising model: 



(ab) a 



(A16) 



The operator has the physical interpretation of creating an 
Ising vortex (vison) on plaquette a. In the deconfined phase, 
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when K/h > x c , the visons are gapped excitations in the 
"paramagnetic" phase of this dual Ising model. We will see 
that the dimerization is related to virtual vison excitations. 

To see this, we obtain a continuum limit of Eq. ( |A16| l, valid 
in the deconfined phase, as follows (qualitatively identical re- 
sults can be obtained in many other ways, for instance by an 
expansion about mean field theory, or by strong coupling ex- 
pansions). It is convenient to work in a path integral formu- 
lation in the basis, and "soften" the spins — > cp a . The 
Euclidean action in the time continuum limit is then 



(ab) 



(A17) 



horizontal bonds, 

D x = {—l) x ' (Si • Si+x — Si+ X ■ Si+2x) (A22) 

~ ( — 1) * { a i,i+x ~ a i+x,i+2x) 

~ (—1) { T a T a+y + T a+x T a+x+y) 

~ (c$i + s<J> 2 )(c$i - s$ 2 ) - + c$ 2 )(s$i - c$ 2 ) 



where in the penultimate line of Eq. ( |A22| i, c = cos7r/8 and 
s = sin7r/8. By a similar calculation, one finds that the verti- 
cal bond dimerization is 

D y = (-1)« (Si • S. l+V - S t+y ■ S i+2 y) ~ 2$!$ 2 . (A23) 

From this we obtain the result 



where k, r and u are phenomenological parameters. In the 
deconfined phase, the fluctuations of ip a are small, and it is 
sufficient to truncate the action to quadratic order. The domi- 
nant fluctuations are those near the minimum of the quadratic 
form. To find them, we must choose a gauge for the frustrated 
dual exchange. It is convenient to make the following choice: 



* = D X + iDy ~ ($1 + i$ 2 



(A24) 



fJ-a,a+y — (~ I)' 1 



f^a,a-\-x 



f. = 1. 



(A18) 



Here we have taken the dual lattice sites to have integer coor- 
dinates. The unit cell in this gauge contains two sites. There- 
fore, Fourier transforming to go to the Bloch basis, we obtain 
the inverse Green's function describing the virtual fluctuations 
of the visons, 



G~ 



:{Kul + r)\-lh ( cos ^ 

v ™ ' \ COS K. 



cos k y cos k x ^ 

- COS ky J 



(A19) 



Here the "magnetic" Brillouin zone is \k x \ < 7r/2, \ k y \ < n. 
The dominant fluctuations, corresponding to the minimum 
eigenvalue of G _1 (= r — A\2K), occur at the two in- 
equivalent values (k x ,k y ) = (0,0) and (k x ,k y ) = (0,7r). 
The corresponding eigenvectors are <jP-' = (cos |,sin |) at 

k = (0, 0) and <j)^ = (sin |, cos |) at k = (0, 7r). Focusing 
on these lowest energy excitations, we therefore write 

Va ~ 0W$ 1 (x a ,y o ) + ^ 2 )(-l)^<i> 2 ( a;a ,y a ),(A2O) 

(i) 

where cpa takes the two values of eigenvector i given above 
when a is on the two distinct sublattices, and $>i(x, y) is a 
slowly-varying continuum field. The bulk effective action is 
then 

S = \ £ / dTdxd y {( d r®if + ^(V*,) 2 + mH%)K21) 

A — 1 o " 



This action describes two degenerate minimum energy vison 
states. It was discussed first to our knowledge in Ref. 46 , in the 
context of frustrated Ising models. It is instructive to express 
the VBS order parameter in terms of If we consider the 



The gauge invariant combination on the right hand side can 
thus be identified as the familiar complex VBS order param- 
eter This result, and the action Eq. ( |A21| >, have been 
obtained many times for quantum spin-1/2 systems on the 
square lattice. Indeed, both are largely independent of the mi- 
croscopic model, and give the minimal set of excitations and 
their properties gives only the assumptions of Z 2 topological 
order in the ground state and half-integer spin per unit cell. 
It would be interesting to understand if other dimer patterns 
could in principle arise, if the low energy vison states were 
selected from a different projective symmetry group. 47 In two 
dimensions, in the Z 2 QSL phase, there is no VBS order, so 
the visons are gapped and the VBS order parameter also is 
uncondensed, correspondingly. 

We now consider the finite-size effects. Taking periodic 
boundary conditions on ip a in the y direction imposes, us- 
ing Eq. (A20i, periodic boundary conditions on <t>i but anti- 
periodic boundary conditions on $ 2 when L y is odd. The 
latter result can be readily understood in terms of the VBS or- 
der parameter: on an odd-leg cylinder, the vertical component 
D y is frustrated (staggering of rows of dimers does not "fit") 
and should be antiperiodic, which requires — > vj/* under the 
circuit around the cylinder, consistent with the anti-periodic 
boundary conditions on $ 2 . Since the visons are gapped, the 
antiperiodic boundary condition gives an exponentially small 
effect in the thermodynamic limit, but it is non-zero and can 
be readily calculated. 

Regardless of boundary conditions, because $i and <i> 2 are 
decoupled in Eq. ( |A21| i, (D y ) = 0, so there is no VBS order 
of the vertical bonds. The horizontal component, however, is 
non-zero when L y is odd, so that the fields $i and $ 2 are 
slightly inequivalent due to the boundary conditions: 



(A25) 



duj n dk x 
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ui% + v 2 k 2 + m 2 



E 



v 2 k 2 
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where the first sum is over "periodic" momenta k y — 
2irn/L y , and the second sum (with the prime) is over "an- 
tiperiodic" momenta k y = 2n(n + 1/2) /L y , with integer n. 
To proceed, we first perform the frequency integration and 
then use the Poisson resummation formula to obtain 



(D, 



2 v ^ f dk x f 
J 2tT J 



dk y cos[(2p + l)k y Ly] 
2n \Jv 2 k 2 + to 2 



(A26) 



dk, 



TTKV z — ' / 2lT 
p=0 



E 



M(k x )K [(2p+l)M(k x )L y /v) 



where we carried out the k y integration in the last line, and 
defined M(k x ) = yjm? + v 2 k x . For large L y , the asymp- 
totic form of the Bessel function can be used, Ko(z) ~ 
y/n/2ze~~ z , and the dimerization is dominated by the p = 
term and the region vk x <C m : 



(D-) —J^— I — me- mL y /v e- vk * L y /2m 



2 j ttv f dk x 

TKV ' 

V2m 



ttkv y mLy J 2tt 



itkvL 



-mLy/v 



(All) 



As promised, we obtain exponential decay of the dimeriza- 
tion, and in this case a prediction for the prefactor. The 
physics of this derivation is transparent: virtual fluctuations 
of Z2 vortices which propagate about the cylinder lead di- 
rectly to the dimerization. In this way we immediately see 
that this effect is universal for Z 2 QSLs on the square lattice 
with S — 1/2 spins. 

Let us conclude this subsection with one remark on the 
dimer correlation lengths. The static dimerization on cylin- 
ders with odd circumference decays with an apparent corre- 
lation length £ = v/m. This is not the same length which 
appears in the dimer-dimer correlation function. The latter is 
obtained from correlation functions of fy, given in Eq. ( |A24| i. 
Because the dimer order parameter \& is quadratic in the vison 
fields <£>i, and the <£>i are Gaussian distributed, by Wick's theo- 
rem the dimer-dimer correlation functions are squares of vison 
Green's functions. Consequently, the exponential decay of the 
dimer-dimer correlation function, which defines the standard 
dimer correlation length £, is twice as fast, i.e. 



(A28) 



This behavior has indeed been observed in the numerical stud- 
ies in the main text. 



2. Ground state degeneracy 

It is well-known that the Z2 spin liquid has degenerate 
ground states in the thermodynamic limit on a cylinder or 
torus. For the cylindrical geometry studied here, two states are 
expected. Here we would like to understand the scaling of the 
gap between these two states, and also better understand their 
character. We will see that, as discussed e.g. in RefF^J, that 
the presence of gapped spin excitations (which carry non-zero 



electric gauge charge) makes a qualitative difference in these 
properties. This means that models neglecting these excita- 
tions, in particular the very popular quantum dimer models, 
actually give incorrect or non-generic scaling for the finite- 
size quasi-degenerate gap. 

Consider first the pure gauge theory, Eq. ( |A3[ ), in which 
coupling to matter fields is neglected. The ground state degree 
of freedom may be regarded as the presence or absence of a 
vison through the hole in the cylinder. The presence of the 
vison itself is measured by the Wilson loop operator around 
the cylinder, 



W 



JT a xy;xy+l- 
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A state with a Z2 vortex in it has W = — 1 and without has 
W = 1. However, the ground state will not be an eigenstate 
of W. In fact, consider the conjugate operator 



n 



xy;xy+l- 



(A30) 



This operator commutes with H defined in Eq. ([A3), and so 
is a constant of the motion. The two degenerate ground states 
have Q — ±1 (we can pick any y, since others are related by 
Eq. ( |Al) ). Note that WQ = -QW, so an eigenstate of Q 
is a symmetric or antisymmetric combination of the W (vi- 
son) eigenstates. This indicates physically that the vison may 
tunnel through the cylinder, by moving (virtually) through the 
entire long length L y from one end to another, thereby con- 
necting the W = 1 and W = — 1 states. The tunneling am- 
plitude for this process is naturally expected to be exponential 
in the length of the event, so we postulate that the gap in this 
case is t v ~ e - L */t* . This has been shown explicitly in many 
places in the literature. 

This result is generic for the pure Z2 gauge theory, and con- 
tinues to hold even if longer (but finite) range plaquette and 
electric field terms are included. It relies only on the fact that 
Q does not create any physical gauge flux through finite pla- 
quettes. However, if a matter field (i.e. the spinons) is present, 
the result is modified. To see this, let us imagine more care- 
fully integrating out the spinons in going from Eq. ( |A1[ ) to 
Eq. |A3), for the case of a cylinder of finite circumference. 
Then we will obtain not only contributions from small loops 
(which renormalize K etc.), but also, occuring first at 0(t Ly ), 
contributions from loops which encircle the cylinder. Keeping 
just the leading of these terms, we have the slight modification 
of Eq. ([A3) 
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(A31) 

where we expect t s ~ e~ Ly ^ y , which physically is related 
to the amplitude for a virtual spinon to encircle the cylinder. 
Note that in this case Q no longer commutes with H, and the 
nature of the eigenstates is no longer clear. Now if we assume 
t s <C K, h and that for t s — we are in the deconfined Z2 
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phase, we can project the Hamiltonian on the low-energy sec- 
tor of the pure gauge theory, i.e. the two level system of the 
quasi-degenerate states. Then we obtain the effective Hamil- 
tonian, written in a pseudo-spin notation in which /j, z = =pl 
correspond to the vison/no-vison states: 

H deg = -t vf i x ~t s L x ^ z . (A32) 

Since t v ~ e -W<5* an( j t s ^ e - L v/£y^ the nature of the 
ground state depends crucially on the aspect ratio of the cylin- 
der. For a fat cylinder, with small L x /L y , for which t s <C t v , 
the eigenstates will be like those of the pure gauge theory, and 
the gap will be exponentially small in L x . 

However, for a "long" cylinder, with larger L x /L y , the 
gap will be exponential instead in L y . Indeed, strictly in the 
limit of large L x and L y fixed, the higher energy state can 
no longer be regarded as quasi-degenerate: its energy, rela- 
tive to the ground state, grows linearly with L x , and so other 
states with local, non-topological excitations will have lower 
energy. The conditions for t s to dominate are much less re- 
strictive than this, however, requiring only t s L x ^> t v , or 
exp(L x /£ x — Ly/^y) ^> 1/L X . In this limit, the ground state 
is an approximate eigenstate of /i z , i.e. a state of definite vison 



number. Because of the quasi-one-dimensional nature of the 
DMRG technique, in the most effective regime of this tech- 
nique, this is the expected form of the ground state. Note 
again that this regime is missed by the pure gauge theory and 
also the quantum dimer model. 

The nature of the absolute ground state obtained by DMRG 
has implications for the entanglement entropy. As shown re- 
cently by Zhang et aF^I the topological entanglement entropy 
for a cut with non-trivial topology actually depends upon the 
choice of quasi-degenerate wavefunction. The cylindrical cut 
studied here is precisely such a cut. The results of RefP^I im- 
ply that the topological entanglement entropy reaches its max- 
imum and universal value (of — In 2) when the ground state is 
a vison eigenstate, and takes a smaller (in magitude) value for 
other superpositions of states, vanishing for the case of a vison 
superposition, as is obtained in the absence of spinons. Thus 
the result of our numerical study in the main text, in which 
we found rough agreement with the — In 2 value for the topo- 
logical entanglement entropy, in fact is evidence for such a 
vison eigenstate in the numerics, consistent with the predicted 
effects of virtual spin fluctuations. 



